% This file provides a Monte Carlo exercise to study the suitability of the
% SLP methodology. It generates Figure D1, presented in Appendix D in the
% paper. 

clear; close all; clc;

%% Generate the True DGP. 

A           =   [0.7 0.1;
                0.2 0.8];
B           =   [1 0;0.2 1];

H           =   9;                      % Number of horizons.
IRF         =   zeros(2,H);
for jx = 1:H
    IRF(:,jx)   =   A^(jx-1)*B(:,1);
end

%% Generate Fake Data.  

N           =   1000;                   % Number of Monte Carlo Simulations.
TT          =   200;                    % Number of periods.
burn        =   100;                    % Number of periods to burn-in. 
T           =   TT + burn;              % Total number of periods to generate.

% Options for S/LP
p           =   1;                      % Number of lags.
lambda      =   100;                    % Value of penalty.
r           =   3;                      % (r-1) = order of the limit polynomial 
                                        % (so r=3 implies the IR is shrunk 
                                        % towards a quadratic polynomial).

randn('seed',727);
e           =   1*randn(2,T,N);
Zsim        =   zeros(2,T,N);
IRFmc       =   zeros(2,H,N);

for jx = 1:N
    Zsim(:,1,jx)    =   B*e(:,1,jx);
    for ix = 2:T
        Zsim(:,ix,jx)    =   A*Zsim(:,ix-1,jx) + B*e(:,ix,jx);
    end
    [AA,SIGMA,U,V]  =   olsvarc(Zsim(:,burn+1:end,jx)',p);
    B2              =   chol(SIGMA(1:2,1:2));
    B2              =   B2';
    [IRFM]          =   IRFVAR(AA,B2,p,H-1);
    IRFmc(:,:,jx)   =   IRFM(1:2,:);
end

% Compute average across runs:

IRFmca      =   zeros(2,H);

for kx = 1:2
    for ix = 1:H
        IRFmca(kx,ix)   =   mean(IRFmc(kx,ix,:));
    end
end

irflp       =   zeros(N,H);
irfslp      =   zeros(N,H);

% Now do a local projection of Y on X:
for jx = 1:N
    x           =   Zsim(1,burn+1+p:end,jx)';
    y           =   Zsim(2,burn+1+p:end,jx)';
    w           =   [Zsim(1,burn+1:end-1,jx)' Zsim(2,burn+1:end-1,jx)'];
    h1          =   1;
    lp          =   locproj_var(y,x,w,h1,H,'reg'); 
    irflp(jx,:) =   lp.IR(2:H+1)'; 
    % Extract the IRF from local projection: 
    slp         =   locproj_var(y,x,w,h1,H,'smooth',r,lambda); 
    irfslp(jx,:) =  slp.IR(2:H+1)';
end

irflpm      =   zeros(1,H);
irfslpm     =   zeros(1,H);
for jx = 1:H
    irflpm(1,jx) =  mean(irflp(:,jx));
    irfslpm(1,jx)=  mean(irfslp(:,jx));
end

%% Creating figure.

t = 0:H-1;

figure
subplot(2,2,1)
plot(t,IRF(2,:),'-k',t,IRFmca(2,:),'--k',t,irflpm,':k',t,irfslpm,'-b','Linewidth',1.5)
title('Average','Interpreter','latex','fontSize',14)
h = legend('True','VAR','LP','SLP','Location','SouthEast');
set(h,'Interpreter','latex');
xlabel('Horizon','Interpreter','latex','fontSize',12);
grid on
subplot(2,2,2)
plot(t,IRF(2,:),'-k',t,IRFmc(2,:,20),'--k',t,irflp(20,:),':k',t,irfslp(20,:),'-b','Linewidth',1.5)
title('Draw 20','Interpreter','latex','fontSize',14);
xlabel('Horizon','Interpreter','latex','fontSize',12);
grid on
subplot(2,2,3)
plot(t,IRF(2,:),'-k',t,IRFmc(2,:,200),'--k',t,irflp(200,:),':k',t,irfslp(200,:),'-b','Linewidth',1.5)
title('Draw 200','Interpreter','latex','fontSize',14);
xlabel('Horizon','Interpreter','latex','fontSize',12);
grid on
subplot(2,2,4)
plot(t,IRF(2,:),'-k',t,IRFmc(2,:,650),'--k',t,irflp(650,:),':k',t,irfslp(650,:),'-b','Linewidth',1.5)
title('Draw 650','Interpreter','latex','fontSize',14);
xlabel('Horizon','Interpreter','latex','fontSize',12);
grid on



